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Abstract 

Nontrivial benchmark solutions are developed for 
the galactic ion transport (GIT) equations in the 
straight-ahead approximation. These equations are 
used to predict potential radiation hazards in the up- 
per atmosphere and in space. Two levels of difficulty 
are considered: (1) energy independent and (2) spa- 
tially independent. 

The analysis emphasizes analytical methods never 
before applied to the GIT equations. Most of the 
representations derived have been numerically imple- 
mented and compared with more approximate calcu- 
lations. Accurate ion fluxes (to 3 to 5 digits) are 
obtained for nontrivial sources. For monoenergetic 
beams, both accurate doses and fluxes are found. 
The benchmarks presented herein are useful in as- 
sessing the accuracy of transport algorithms designed 
to accommodate more complex radiation protection 
problems. In addition, these solutions can provide 
fast and accurate assessments of relatively simple 
shield configurations. 

1.0. Introduction 

As mankind turns its attention toward establish- 
ing orbiting space stations to serve as bases for fu- 
ture space exploration, the protection of personnel 
against radiation has become a relevant issue. In the 
upper atmosphere and in space, high-energy heavy 
ions originating in deep space or in the Sun are a 
major source of penetrating particles. Not only is 
the unscattered contribution of concern, but also the 
secondary particles generated from direct nuclear or 
coulombic collisions can account for a significant frac- 
tion of the absorbed dose. To ensure proper shielding 
of space bases and interplanetary vehicles, predic- 
tive dose and flux computations must be developed. 
Monte Carlo methods that numerically simulate par- 
ticle motion are important calculational tools; un- 
fortunately, however, these methods require exten- 
sive computational resources for adequate accuracy 
and, therefore, cannot be routinely applied. With 
recent advances in numerical methods and compu- 
tational strategies, deterministic methods character- 
ized by the Boltzmann equation have become com- 
petitive to Monte Carlo simulations. For this reason, 
an effort has been initiated within NASA to develop 
reliable, multidimensional, deterministic methods for 
shield design. 

One component of this effort has been directed to- 
ward the assessment of the accuracy of proposed de- 
terministic algorithms. The “partial” verification of a 
given algorithm can be obtained by comparison with 
standard solutions of the governing equations. These 
solutions, referred to as benchmarks , are highly ac- 


curate numerical evaluations of simplified transport 
problems which, nevertheless, contain the physical 
features of the basic transport processes. The mo- 
tivation behind such comparisons is that algorithms 
developed for realistic problems must also yield reli- 
able results for these simple problems as well. Com- 
parisons with benchmark solutions are just one com- 
ponent in the verification process. Others include 
particle conservation, following proper intuitive phys- 
ical trends, and comparisons of results with those ob- 
tained from alternative algorithms designed for sim- 
ilar transport problems. Even with these additional 
components, only partial verification is achieved be- 
cause each application is unique. 

In this report, benchmarks for the straight-ahead 
Boltzmann equation describing ion transport in mat- 
ter are developed. Particular emphasis is placed on 
the analytical techniques used to solve the coupled 
ordinary differential transport equations. The ap- 
proach taken is to investigate simplified problems, 
eventually leading to consideration of the full set of 
galactic ion transport (GIT) equations. The current 
work is devoted to two problems, energy-independent 
and spatially independent problems. A full analy- 
sis is performed for each problem, with results for 
physically interesting cases included. In this way, 
the knowledge gained at each stage of the investi- 
gation can be employed in the next. Although the 
major emphasis herein is on the development of ac- 
curate numerical evaluations, the solutions obtained 
can also be used to provide an approximate charac- 
terization for realistic transport situations. 

2.0 1 . Galactic Ion Transport (GIT) 

Equations 

Because of the high energy of the galactic ions 
(10 to 10 7 MeV), the straight-ahead approximation 
can be introduced into the Boltzmann equation with 
a high degree of confidence. In this approximation, 
ions are not angularly deflected; and, as the collid- 
ing ions break up in nuclear fragmentation, the frag- 
ments continue in the incident ion direction. Thus, 
for ions of charge number j, the appropriate trans- 
port equation is (ref. 1) 


' d_ 
dx 


d 

dE 


Sj(E) + 


4>j{x,E) 


J 

k=j + 1 


(la) 


where <f>j is the flux of the jth ion at position x with 
energy per nucleon of E. (A list of symbols appears 
after the references.) The macroscopic absorption 
cross section is mjk is the multiplicity of ion 
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j produced in collision with ion k , and Sj is the 
absolute value of the change in energy E per unit 
distance traveled (i.e., the stopping power). The 
added assumption that target ion fragmentation can 
be neglected (ref. 1) is made in equation (la). For 
most (but not all) cases studied, a homogeneous 
semiinfinite medium with a source at the free surface, 

4> j {0,E) = f j {E) (lb) 


(^ +<7j) M x ) = o (6) 

and y v 

MO) = 1 (7) 

The appropriate solution to equations (6) and (7) is 
therefore 

<j)j{x) = exp(-o-jx) (8) 

For j = J — 1 we have 


is assumed. 

As each benchmark is developed, additional as- 
sumptions appropriate to that benchmark are made. 
Also, in certain instances the restriction to a homo- 
geneous shield is relaxed. 

3,0. Energy-Independent Benchmark 

3.1. Theory 

We may obtain the energy- independent bench- 
mark from the Boltzmann equations (la) and (lb) 
by first assuming that the cross sections and frag- 
mentation multiplicities are constants (independent 
of E). Equations (1) are then integrated over all en- 
ergies to yield 



4>j{ 0) = 6jj (3) 


(^ + <7/-l) 4>j-i( x ) = 0J-l,J<l>j( x ) (9) 

with 

<t>j-\(o) = h-\,j = o (io) 

Now equation (9) is of the general form 


JL+P(x)y = Q(: v) 


( 11 ) 


which can be solved using an integrating factor as 
y = exp [-/ P{x) dx ]{/ <5(ar)exp U PMd * M 

+ Cexp^— J P(a;)dx| 


( 12 ) 


Because 


and 


y = <t>j - iW 

P(x) = CT/.! 


where the energy- independent flux is given by 


Q( x ) = Pj-i,jM x ) = Pj-i,j exp(— < tjx) 


and 


4>{x) 


-r 

JO 


4>j(x , E) dE 


fijk ~~ ^jk^k 


(4) 

(5) 


A source of only J ions on the free surface has 
been specified. Equations (2) and (3) are solved 
with three different methods: (1) direct analytical 
solution, (2) Laplace transformation, and (3) power 
series expansion. Equation (2) is also applicable to 
a heterogeneous medium with spatially varying cross 
sections. This is accomplished by interpreting x as 
an optical thickness rather than as a linear thickness. 


3.1.1. Direct analytical solution 

Since equation (2) is nothing more than a set of 
coupled, ordinary differential equations, the general 
solution is a linear combination of exponentials. For 
example, for j — J we have equations (2) and (3) 
reducing to 


the solution of equation (9) is 
4>j- 1 (*) 

= exp (~<tj-ix) ^ J Pj-i,jM x ) exp {<tj-\x) dx 

+ C exp(-«7 t /_ 1 x) 

= exp(— <Tj_ 1 rr)/3 < /_ 1 j J 

x exp(-crjx) exp (<Jj~ix) dx 
+ C exp(-arj_ix) 

= exp (-<Jj-ix)f3j- h j J exp(<r t /_ 1 — a j)xdx 
+ C exp(-aj_!x) 

= exp{-aj_ix) exp [{crj^ - < tj)x ) 

a J - 1 -°J 

+ C exp(-aj^ix) 


°j~\ - °j 


exp(— ojx) + C exp(— aj^^x) (13) 
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Now, applying equation (10) to (13) yields 


to equations (2) and (3) yields 


4>j- 1(0) = o = 


— — +c 

a J - 1 - a J 


hip) = 


i_ 

p + CTj 


s jj + 5Z fljkfikip) 

k=j + 1 


( 19 ) 


which implies that 

c = 

aj_i - gj 


Equation (13) thus becomes simply 

<t>j-i{x) = /3j ~ — [exp(-o-j-x) - exp(-ffj_ 1 x)] 

^j-i - 

(14) 

From a constructive inductive proof, the general 
solution to equation (2) is 

l 

i = 0 

(0</<J-l) (15) 


where the coefficients are given recursively by 


a 0,J = 1 

a i,J-l = 2 ~ “ a i,J-k' 


(16a) 


/-i 

a l, J-l = ~ 

i = 0 


(0 < i < l - 1) (16b) 
(16c) 


The last expression for is required to satisfy 

the boundary condition (eq. (3)) 


4>j( 0) = 0 (l<j<J-l) (17) 


Because of the recurrence relation given by equa- 
tions (16), the expression for cf>j-i(x) in equation (15) 
is susceptible to round-off error when J is large, and 
this error could be a significant numerical limitation. 
This direct analytical method is also applied to the 
energy-dependent benchmark. 


Thus, an algebraic recurrence relation for the trans- 
formed flux has been established. This recurrence 
relation is most conveniently used in a numerical in- 
version procedure based on a decomposition of the 
Bromwich inversion integral into an infinite series, 
with application of the Euler-Knapp transformation 
to accelerate its convergence. When J is large, how- 
ever, this procedure is expected to require excessive 
computational time as a result of the large number 
of transformed function evaluations ((f) j(p)) required 
for the numerical integration. 


3.1.3. Power series expansion 


The least complicated solution method is based 
on a simple power (Taylor) series expansion of the 
form 


= exp(-rx) ~ 7 h n,j 

n=0 U ' 


(20) 


Introducing this representation into equation (2) 
yields 


d(f) 


— x n 


-at = -re*p(-r*)E-|Jvi 


77=0 

x n- 1 


+ exp(-Tx) J2 (~ _ iy h n,j 


so that 

°° „n- 1 




x" * _ ^ x 

^ (n-l )! n,j “ nf 

71=0 V / 77=0 


(gj T)h n j 


+ 0jk^n,k 

k=j + 1 


Equating coefficients of like powers of x yields the 
recurrence relation 


J 

^n+l J = ~~ T)hn,j “k ^ 0jk^"n,k (^1) 

k=j+l 


3.1.2. Numerical Laplace transform inversion 


From equation (20) we note that 


Applying the Laplace transform 


W) = djj = h o ,j ( 22 ) 


roc 

<f>j(p) = dx exp (—px)cf)j(x) 

J 0 


(18) 


The introduction of the multiplicative exponential 
factor exp(— Tx) in equation (20) is used to speed up 
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convergence of the power series in neutron transport 
calculations, but is not required here (T = 0). 

A variant of the above expansion which can be 
used to furt her speed up convergence for large x is ob- 
tained by partitioning the slab into intervals. These 
intervals may be chosen to correspond to different 
shield materials or chosen for computational conve- 
nience. For the ith interval, defined as < x < a;, 
the following expansion is made: 

= t (23) 

n=0 

The recurrence relations for h l n found in an anal- 
ogous manner to those obtained in equations (21) 
and (22). are 

4j = *5 -1 («i-i) (24) 

and 

J 

K + i = £ Wuk (25) 

k=j+l 

Equation (24) merely asserts that the flux incident 
on the ith slab is that which exits the i — 1 slab. 
Again, one must be alert to potential round-off error 
when using equations (20) to (25). 

3.2. Numerical Implementation 

In order to numerically compare the above solu- 
tion methods, the following simplified nuclear model 
is assumed: 


■j - V 0 j 213 


(26) 

2 

(k > j) 1 

(27) 

k- 1 

0 

(* < J) i 



The choice of aj is based on the liquid drop model of 
nuclear physics. The multiplicities are chosen so as 
to conserve charge in each interaction. 

3.2.1. Solution comparisons 

Table I displays a comparison of the total flux 

J 

«M*) = 5>j(x) ( 28 ) 

j = 1 

calculated using analytical and power series expan- 
sion methods for incident manganese ions (J — 25) 


colliding with an air shield ( <j 0 — 0.01247 cm 2 /g) 
100 g/cm 2 thick. 

In the power series method, adequate convergence 
is obtained when the sum of three individual terms 
yields a relative error less than some desired amount 
In the table, the total fluxes for four relative er- 
rors are shown compared with the fluxes for the an- 
alytical solution. The underlined digits are those in 
disagreement with the analytical solution. All power 
series calculations are within the specified relative 
errors. The power series method, however, takes ap- 
proximately twice the computational time of the ana- 
lytical method. In most respects the two methods are 
the same, except that the power series is more easily 
applied to heterogeneous shields. Also, and possibly 
more importantly, if round-off error becomes a prob- 
lem, there is a readily available remedy for the power 
series expansion (multilayer treatment) but not for 
the analytical solution method. 

A comparison of the analytical solution with 
the numerical Laplace transform inversion (at x = 
100 g/cm 2 ) for different relative errors ep is displayed 
in table II. For 5-digit accuracy, the numerical inver- 
sion approach is inferior to the analytical or power 
series solution methods as it requires about 800 times 
the computational effort. Because of this computa- 
tional inefficiency, the numerical inversion is not con- 
sidered further for this benchmark. 

Several physical trends are displayed in fig- 
ures 1(a) and 1(b). Figure 1(a) shows the variation 
with x of the 25 species. As anticipated, the flux in- 
creases with decreasing j . This, of course, is a result 
of the increased fragmentation with increasing depth 
x as the ions penetrate farther and undergo more col- 
lisions. In figure 1(b), the variation of 4>t{ x ) with x 
is shown. Beyond x — 100 g/cm 2 , ftp begins to de- 
crease and eventually becomes zero at infinity (large 
x). This is a physical requirement since all species 
have nonzero absorption cross sections, a fact which 
implies that each species must eventually disappear. 

3.2.2. Multilayer formulation 

The flux at 50 locations within a 100 g/cm 2 atmo- 
spheric shield was determined for both a single layer 
and for 5 layers each with 10 locations. The compu- 
tational time for the five-layer calculation was about 
three-tenths of the time for the single layer, which in- 
dicates a clear advantage to using the multilayer for- 
mulation. As displayed in figures 2(a) and 2(b), the 
major advantage of the multilayer formulation is the 
more rapid convergence in each new layer. In these 
figures, the number of terms required for conver- 
gence N conv (for ep — 10~ 4 ) for the flux of the last- 
generated ion </>i (x) is shown for both single-layer and 
five-layer configurations. A maximum of 10 terms are 
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required for the 5-layer calculation, whereas 23 terms 
are required for the single- layer calculation. This 
accounts for the previously mentioned reduction in 
computational time. This time advantage will, how- 
ever, become a time disadvantage as the number of 
layers increases because of the associated increase in 
overhead. Another, more subtle, advantage of the 
multilayered formulation is the mitigation of round- 
off error by the accelerated convergence. If round-off 
error is suspected, or if the series does not converge 
within a specified number of terms, then more layers 
can be introduced to reduce the number of terms per 
layer required for convergence. Thus, fewer iterations 
of the recurrence relations are needed to provide more 
accurate values of h n j to be used in the summation. 


3.2.3. Comparison with other methods 

A natural numerical approach to solution of equa- 
tions (2) and (3) is the standard implicit numerical 
finite-difference scheme. If the shield is partitioned 
into N uniform intervals of width Ax and the deriva- 
tive in equation (2) is approximated at the nth inter- 
val boundary by 


d^jjx) _ - 4> r j 1 

dx Ax 


(29) 


then the following recurrence relation ensues: 

<t>) = Sjj (30) 


*> + ' = rr+ ($ + A +/^ +I ) (3D 

In table III, the total flux (eq. (28)) obtained from 
this implicit discretized formulation is compared to 
the analytical result for decreasing Ax (increasing 
N). Surprisingly, better than 5-percent accuracy is 
obtained for 20 intervals and 0.1-percent accuracy for 
160 intervals. In addition, the correct result is ap- 
proached as Ax decreases. Thus, a finite-difference 
formulation can provide an acceptable numerical so- 
lution to the energy-independent GIT equations. 

A second comparison is now made with a mul- 
tiple collision (or perturbation) formulation devel- 
oped in reference 1. The estimated percent differ- 
ences between the multiple collision and analytical 
solutions are listed in table IV. Based upon the last 
neglected term, it is estimated that the multiple col- 
lision formulation is accurate to within 10 percent. 
This is clearly verified by the analytical benchmark 
comparison. 


3.2.4 * Selected case studies 

3.2.4- 1- Generation of ion J + 1. Because of 
proton capture reactions, there is a finite probability 
that ion J + 1 can be formed from ion J. To account 
for this, an equation for ion J+l must be considered. 
Therefore, equation (2) becomes 

+ aj+1 ) 4>J+x{x) = c+<t>j{x) (32) 

and 



where the sum in equation (33) now includes J + 1, 
since it can fragment into lighter ions (I < j < J) 
after it is formed. Equations (32) and (33) are most 
easily solved with the power series expansion to give 



h 0,j = <t>) 

(34a) 

K,j + 1 = 

- ~ a J+ 1K1-1J+1 + C+Ki-i,j 

(34b) 

Kj = 

-ajh l n _ X J + 0 jj+\h l n -\j + i 

(34c) 


J+l 


^ n,j = 

+ E fykKi-i,k 

(34d) 


k=j+ 1 



Figures 3 and 4 show the influence of C+ on the total 
flux equation (28) with J replaced by J + 1 and on 
</>j + l(x). The most probable value for C+ is about 
0.01 so that the generation of ion J T 1 does not 
have a significant impact on f>p. If C+ increases to 
about 0.1, however, then the production of ion J + l 
can become important. Figure 4 shows the expected 
reduction in <t>j+\ profile with decreasing C+. 

3.2. 4-2. Multiple shields. The total flux for a 
composite shield of 10 layers with differing density 
factors relative to air are displayed in figure 5. As 
is apparent in the low-density shields (5, 7, and 10), 
the total flux remains nearly constant because of the 
reduced collision probability. For the high-density 
shields (4, 6, and 8), the flux changes rapidly be- 
cause of the increased fragmentation probability. In 
shield 4, the flux initially increases (because of frag- 
mentation) but then decreases because of the pre- 
dominance of absorption. 

4.0. Spatially Independent Benchmark 

4.1. Theory 

A spatially independent benchmark is easily de- 
rived by considering cross sections and multiplicities 
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as spatially constant and treating the boundary con- 
dition as a volume source. For a monoenergetic 
source of ions J, equations (la) and (lb) become 



4>j(x,E) 


J 



= 'Z Pjkity&kiX’E) 



fc=j + l 

+ 6(x)6(E)6jj 


(35a) 

<t>j(0,E) = 0 


(35b) 


Integrating over the half-space (0 — ► oc) and defining 

roc 

cf)j(E) = J cf)j(x,E)dx (36) 

yields the desired spatially independent ion transport 
equations 


and 

4>j{ 0) = bjL (40b) 

vj 

where the path length is given by 

r.Eo dE f 

s{e ' e ° )= L wm (41) 

which, from the definition of the flux distribution, 
implies 

<f)j{x,s) = S p (E)4>j(x, E) (42) 

The major difficulty in applying the above transfor- 
mations is the use of equation (41) to determine the 
path length s from the energy E since an integral 
must be inverted to return to the energy-dependent 
flux of equation (42). This added inconvenience, 
however, is readily acceptable because of the enor- 
mous simplification offered by the energy and path 
length transformation. Indeed, equation (40) now 
closely resembles the energy- independent case, and 
similar solution methods can therefore be adopted. 


-±S,(E) + a,(E) 


4>,{E)= 

k = j + 


0 jk (E)<t> k (E) (37a) 


and 

^ = sjik (37b) 

The volume source has been transformed into an 
initial condition in energy by our integrating over an 
infinitesimal interval about E 0 and noting that </>(E) 
at the upper bound is zero since upscattering is not 
allowed. 

At this point, it becomes convenient to introduce 
a well-known scaling law between stopping powers at 
high energies: 


Sj(E) = v'jSp(E) (38) 


where 



(39) 


and S p is the proton stopping power, Sj is the stop- 
ping power of the jth ion, Zj is the charge number 
of the j th ion, and Aj is the atomic mass of the jth 
ion. With this substitution and variables changed 
from energy to proton path length traveled at a given 
energy (s(E,E 0 )), equations (37a) and (37b) become 



+ CTj(s) 


4>j(s) 


j 

0jk(s)4>k(s) 


k=j + 1 


(40a) 


4.1-1- Analytical flux solution and dose 

As before, an analytical solution can be found for 
constant multiplicities and cross sections. By casting 
equations (40a) and (40b) in the form of equations (2) 
and (3), the solution becomes (from section 3.1.1) 


i = 0 

(1 = 0, 1, . . . , J — 1) (43) 


where 


- v, 


(44a) 


a i,J-l = 


[ (aj-i/vj-i) ~ (vj-i/vj-i) 

1 J-l+1 

x 5Z &J-l,k a i,k 


vj-l 


k—J—i 


l - 1 
i = 0 


(0 < i < / — 1) (44b) 
(44c) 


To treat explicit energy dependence of aj and 
mjk, a piecewise constant assumption, or interval 
formulation, similar to the multigroup approximation 
of neutron transport theory is made. The path length 
(or energy) variable is partitioned, with constant 
transport properties assumed within each interval 
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(or group). Equations (40a) and (40b) can then be 
written in each interval as 

("j^ + a f) + 4> r j(s) = E Pjk<t>k ( s ) ( 45a ) 

k=j+l 

and 

(f)j (fl r _i ) = (f>j (fl r _i) ( a r—l < & < a r) (45b) 

where continuity of the flux is maintained across each 
interval interface. The quantities a r - and /3j k are the 
appropriate transport coefficients averaged over the 
rth interval. Again, the solution is of the form given 
in equation (43) as (for the rth interval) 

l 

4> T j-i{s) — E a i,J-l eM-Vj-is/vj-i) 
i = 0 

(/ = 0, 1, 7-1) (46) 

with 


Substituting equation (46) for each interval gives 



x exp [— — a T -\)lvj_i\ (52) 

To perform the indicated integrations, an explicit 
relation between s and E (eq. (41)) must be specified. 
The analysis is somewhat simplified, however, if the 
following analytical range-energy relations (ref. 2) are 
used: 

s{E 0 ,E) = R(E 0 )-R{E) (53) 

where 

R{E) = a 0 \n{l + a l E n °) (54) 

The term R(E 0 ) is actually the maximum proton 
range. Inverting equation (53) for E gives 

E = [(exp { [R(E 0 ) - s] /a 0 } - l) / ai ] 1/n ° (55) 


a 




Gfc) 


J-l+1 

X E fo < * < / - 1) (47) 

k—J—i 

for each interval r. For i = l, we have 

/-I 

a UJ-l = 1) - E a i,J-l ( 48 ) 

i = 0 


This formulation is completely analogous to the spa- 
tial multilayer formulation discussed in section 3.1.3. 

With this particular form of the analytical solu- 
tion, the evaluation of dose becomes relatively sim- 
ple. The usual expression for the energy deposited 
by ion j is 


The expression for Dj given by equation (52) can 
therefore be simplified to 



4-1-2. Power series expansion 

If piecewise constant transport properties are as- 
sumed, a formulation identical to the multilayer so- 
lution developed in section 3.1.3 can be employed. 
Following that same solution procedure, we assume 
the expansion 


f E o 

Dj(E) = Aj / dE SjiEtyiE) (49) 

J 0 

which can also be written as 

[E . 

Dj (E) = vjAj / dE <t>j [s(E, E 0 )\ (50) 

JO 

From the partitioned nature of s, and therefore E , 
we have 



= e ( * :r‘ r ^ < 57 > 

n=0 n - 

which, when substituted into equations (45), yields 
K,j = 4> T j~ l {s r - 1 ) (58) 

and 

K+i = f ( h "j + E /%*<*) (»>) 

Vj V k=j+ 1 J 
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The dose can be determined from equation (51) single region: 


=5. Y hi , /•£,-, 


Dj = "j A J ar-l) n (60) 

n=0 r=l n - ^ 

which, upon insertion of the range-energy relations 
(eqs. (53) and (54)) becomes 

oc N r ur 

= '>*, EE -ft 

x< f E E -- dEi Jl±^lY { 61) 


Not only are these integrals more complicated than 
those for the analytical result (eq. (56)), but also 
there are typically more of them to evaluate, making 
the power series expansion method less attractive for 
calculating doses. 

4.2. Distributed Sources 

For a single ion source of type J at the boundary, 

4 >j (0,E) = 6 jJ f(E) (62) 

The spatially independent benchmark is again ob- 
tained by integrating equation (35a) over x to yield 


- m s j{E) + *j <t>j(E) 


= El PjkhW + SjjfiE) (63) 
k=j + 1 

Changing to the path length variable and using equa- 
tions (38) and (39) we find 


d 1 - J 

'jfa. +<j j( s ) 'M 5 ) = Y Pjkfa(s) + 


where 


fc— j+l 


m = S P (E)f(E) 


With the Green’s function found in the preceding 
two sections, the solution is obtained as the following 
convolution integral: 

M = f du f(s - u)4>j{u) (66) 

and from the analytical and power series solutions, 
one can therefore write the following equations for a 


' 3 rs 

<i>j{s) = Y a ij / duf(s - it) 

i=0 

x exp(— <7j_ju/ vj-i) 
°°h r r$ 

0jOO = Y J T duf{s-u)i 
n] JO 


The integrals in equations (67) and (68) can be costly 
to evaluate numerically; therefore, the method with 
the fewest number of integrals is preferable. In 
most instances, the analytical formulation requires 
the least computational effort. 

Another very effective approximate solution, for 
general distributed sources, can be obtained with the 
multigroup or interval formulation technique. (See 
section 4.1.1.) If the source is averaged over the path 
length interval r(/J), an analytical solution can again 
be found in the form 


j-i{s) - a.j_i + Y a lj~i 
i = 0 

x exp (s-o r _i)] (69) 


where 


0j~i = 


This solution is established by finding the particular 
solution for a constant source and adding it to the 
homogeneous solution. The coefficients and 

the additive term A T j_[ are now more complicated, 
but they can still be found recursively as follows: 

“M = fj- ( 70a ) 


l,J-l - fj-l ~ A J-l - Y 




x fj-l + Y Pj-l,J-i A J- 


with j_ t given by equation (44a) for 0 < i < 
l — 1 within each interval r. This procedure can 
be continued for any polynomial variation of the 
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source, although the algebra becomes increasingly 
more difficult and cumbersome: Because of this 

algebraic complexity, only analytical solutions for 
constant sources are considered. 

The power series expansion, however, is more 
easily extended to other polynomial sources. For 
example, for a quadratic source distribution, 

fj(s) = Q r 0 + Q\(S - a r _i) + Q r 2 (s - a r _!) 2 (71) 
only the term 

(Q r 0 6n, 1 + Qtfn, 2 + QM 

need be added to the right-hand side of equation (59) 
to give the appropriate h r n ■ . 

4.3. Numerical Implementation 

4.3.1. Flux solution comparisons 

Table V provides a comparison of (f>i(E ) as given 
by the analytical and power series expansion solu- 
tions with J — 25 and E 0 = 5 GeV. For the re- 
mainder of this work, the nuclear properties given by 
equations (26) and (27) are utilized unless otherwise 
stated. As shown in the table, the power series solu- 
tion does not uniformly approach the analytical so- 
lution as the truncation error is decreased from 0.01 
to 0.0001. This failure is apparently a result of round- 
off error, as implied by the evaluation of h^j in dou- 
ble precision, for which most of the <j> 1 values compare 
better with the analytical ones. This error is elimi- 
nated when the same energy region is partitioned into 
4 intervals with 10 points per interval. As shown in 
the table, the four-interval double-precision and an- 
alytical results are identical. The reason for this is 
clearly that fewer terms are required for convergence, 
thereby reducing the inaccuracies introduced into the 
calculated coefficients h r n - (section 3.2.2). To reduce 
the prospect of round-off error being a problem in 
determining ajj and h r n ^ these coefficients are cal- 
culated with double-precision arithmetic. 

Figure 6 displays the ion flux variations with en- 
ergy. The flux typically increases with decreasing 
j (for a given energy), as the number of fragmen- 
tation interactions increases the secondary particle 
production. 

4.3.2. Dose calculations 

The dose calculations appear to be less accu- 
rate and reliable than the flux calculations. In ta- 
ble VI, the ion dose factors computed with the ana- 
lytical solution (eq. (56)) are displayed for decreas- 


ing relative error of the numerical Romberg integra- 
tion scheme (ref. 3). Convergence to only two sig- 
nificant figures is achieved. This lack of accuracy, 
resulting from accumulated round-off error of equa- 
tion (56), occurs because subtraction of large num- 
bers eventually yields relatively small numerical val- 
ues. Any variation in these large numbers resulting 
from reductions in the integration error produces dif- 
fering results. For J = 10, however, table VII shows 
that the calculated doses do converge. Therefore, 
we conclude that the dose calculation is useful as a 
tabulated benchmark only if few ions are involved 
(less than 10); otherwise, only graphical comparisons 
should be made. Figure 7 displays such a comparison 
for various initial energies E 0 . As E 0 increases, the 
lighter ions begin to contribute significantly to the 
total dose 

J 

D T = J2 D j < 72 ) 

3 = 1 

since higher E 0 implies that more fragmentations 
occur. 

4.3.3. Distributed sources 

4.3.3. 1. Benchmark source. To illustrate the 
application of the convolution integral (eq. (66)) to 
generating benchmark data, a simple exponential 
source of the form 

f(s) = exp(-a 0 s) (73) 

is considered. In the energy variable, this source 
corresponds to 

f(E) = [S p (£)] _1 Za 0 a 0 (74) 

For this particular source, the integration in equa- 
tion (67) can be analytically performed to yield 

J - j 

Qj(s) - H aij[exp{-Pj_is) 

i = 0 

- exp(— a 0 s)] - a Q ) (75) 

For an exponential source of manganese ions ( J = 25) 
in air, the rapid variation of 4>\ with a Q is displayed 
in figure 8. As shown in figure 9, the energy variation 
of the flux for all species follows that of the source 
ions. Since all calculations involve well-known, on- 
line functions, no convergence problems are experi- 
enced for this source type. Therefore, the convolution 
integral appears to be an efficient way to generate 
reliable results for a class of sources for which the 
integral can be performed analytically. 
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4-3. 3. 2. Solar source. A significant source of 
heavy ions is the Sun. The approximate energy 
dependence is given by (ref. 4) 


/(£) = 


bg 

(6j + E) n 


(76) 


where b 0 , b\, and n are parameters. For both an- 
alytical and power series solution methods, interval 
options are used to evaluate the flux from this source. 
The analytical results for an assumed constant source 
of J = 25 ions within an interval are displayed in 
figures 10(a) and 10(b) for 15 and 40 intervals span- 
ning the energy range of 10 to 10 4 MeV per nucleon. 
There are two points taken per interval. With only 15 
intervals, oscillations in the ion fluxes are clearly no- 
ticeable. When the number of intervals is increased, 
however, the oscillations are reduced (see fig. 10(b)). 
Similar smoothing effects are obtained for the power 
series solutions displayed in figures 11(a) to 11(c). 
Note that the quadratic source representation used 
for figure 11(c) provides fluxes free of oscillations. 

To verify that the interval calculations are indeed 
giving accurate fluxes, the solar source is replaced 
by the exponential source from equation (73). Com- 
parisons between 4>\ calculated with the analytical 
interval methods and those values obtained from the 
benchmark source with a Q = 0 and 10“ 3 are dis- 
played in figures 12(a) and 12(b). As is clearly shown, 
the results for are virtually identical for the two 
solution methods, and this agreement indicates that 
the interval calculation is sufficiently accurate to use 
in estimating fluxes from distributed sources. 

4.3.4- Comparison with a finite- difference 

solution 

A finite-difference scheme similar to the one de- 
veloped for the first benchmark can be constructed. 
This scheme takes the following form for a delta func- 
tion source: 

</>? = 1/Vj (77) 

and 



1 

Vj + As a 9 - 
J J 


v A 1 + As E PjkK 

k=j + 1 


(78) 


and As is determined by dividing the interval 
[0, R(E 0 )] equally: 


j = 1 ion are compared for decreasing As (increas- 
ing N p ). As As decreases, the interval solution does 
indeed approach the analytical solution, with excel- 
lent agreement at high energies (E > 1 GeV/amu). 
For lower energies, the agreement is acceptable but 
not spectacular. Because the majority of the dose is 
contributed by the high-energy flux, the interval cal- 
culation appears to provide an adequate means for 
dose characterization. 

In an attempt to improve on the agreement at 
relatively low energies, a hybrid scheme utilizing 
both the analytical and interval solution methods was 
developed. In this approach, the analytical solution 
is used to obtain the flux for j > Jq and the interval 
solution is used for i <j< Jq — 1. The summation in 
the interval solution (eq. (78)) is split into two sums 

E Pjk^l + E pjkK 

k=j + 1 k=J 0 

where the right-hand summation is obtained from 
the analytical solution. In figure 14, (j)\ is plotted 
for several values of Jo when N p = 10. A signifi- 
cant improvement over the original interval solution 
(Jo = 25) is observed for low energies, even when 
only one analytical flux is used (i.e., Jq = 24). For 
this last case, however, the flux at high energies is less 
accurate than that obtained for the original solution 
(Jq — 25). As Jq approaches 1, the solution continu- 
ally improves. The discrepancy between Jq = 2 and 
Jo = 1 is entirely due to the finite-difference approx- 
imation. Improvement of this can be obtained only 
by increasing N p . 

The complete spectra for all ions using both an- 
alytical and interval solutions ( N p = 100) are dis- 
played in figure 15. These calculations are in very 
good agreement, considering that the finite-difference 
scheme is crude and J is relatively large. 

4.3.5. Energy- dependent absorption cross section 

To demonstrate the energy-dependent absorption 
cross-section option, a functional form, as displayed 
in figure 16, has been assumed. This is similar to 
the actual cross-section shape (ref. 5). All ion cross 
sections are assumed to scale according to 

Oj{E) = j 2 / 3 f c (E) (80) 




R(Eo) 

N p 


(79) 


where N p is the number of intervals. In figures 13(a) 
to 13(d), the interval and analytical fluxes for the 


The fluxes for all ions for cross sections with 
and without the energy dependence are displayed 
in figures 17(a) and 17(b). In general, the high- 
energy portion of the flux (above 1 GeV) is relatively 
unaffected by the energy dependence since the cross 
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sections are essentially constant there. At lower 
energies the fluxes are slightly higher for the energy- 
dependent cross sections, since cfj{E) is larger than 
the corresponding energy-independent values and, 
thus, there is a greater probability of interaction. 

5.0. Concluding Remarks 

Energy-independent and spatially independent 
benchmark solutions to the galactic ion transport 
(GIT) equations have been developed, and accurate 
ion fluxes have been obtained for realistic source rep- 
resentations. Analytical, power series, and Laplace 
transform inversion solution methods were devel- 
oped. The analytical methods are limited only by 
round-off errors. The power series methods are ap- 
proximately half as fast, computationally, as the ana- 
lytical methods, but are more easily altered to correct 
round-off errors. The Laplace transform inversion 
methods are the least useful because of the exces- 
sive computational times necessary for handling the 
large number of different ions in the galactic cosmic 
ray spectrum. For monoenergetic laboratory beams 
both accurate fluxes and doses were found. Examples 


of applications to relatively simple shield configura- 
tions were also presented. The coupled GIT equa- 
tions are partial differential equations which are, by 
their nature, more difficult to solve; nevertheless, an 
analytical solution can be found. 

NASA Langley Research Center 
Hampton, VA 23665-5225 
January 12, 1989 
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7.0. Symbols 

R 

range, g/cm 2 

A j 

mass number of jth ion 

Sj 

stopping power of jth ion, 

ai 

position of zth interval boundary, 


MeV-cm 2 /g 


g/cni 2 

s p 

proton stopping power, MeV-cm 2 /g 

bo ■ b\ 

parameters used in equation (76) 

s(E,E 0 ) 

proton path length, g/cm 2 

C 

arbitrary integration constant 

X 

position, g/cm 2 

C + 

probability of forming ion J + 1 

y 

arbitrary, dependent variable 

Dj 

dose due to jth ion, MeV/g 

a ij 

recursion coefficient 

Df 

total dose, defined in equation (72) 

a 0 ,ai,n 0 

parameters used in equation (54) 

E 

energy per nucleon, MeV /nucleon 

Pjk 

= m jk (J k (eq- (5)) 


or MeV/amu 

r 

convergence factor in equation (20), 

ME) 

energy-dependent absorption cross 
section 

As 

«n 2 /g 

= ^ (eq- (79)) 

Si 

Kij' hruj 
J 

monoenergetic beam source term 
power series coefficients 

Ax 

interval width, g/cm 2 

charge number of incident ion 

bjj 

Kronecker delta 

= jt ( ec l* ( 39 )) 

j 

charge number of jth ion 

lJ j 

m jk 

multiplicity of ion j produced by 
ion k 


macroscopic absorption cross 
section for jth ion, cm 2 /g 

N 

number of intervals 

°o 

cross-section parameter, cm 2 /g 

J^conv 

number of terms for convergence 


flux of jth ion 

P(x).Q(x) 

arbitrary polynomials 

(frj 

Laplace transform of <f)j 

Q r i 

rth interval polynomial source 
coefficient (i = 0, 1,2) 

<Pt 

total flux (eq. (28)) 
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Table I. Total Flux Computations as Function of Atmosphere Thickness 




Total flux (relative to incident flux) calculated with— 





Power series expansion for — 


Thickness, 

g/cm 2 

Analytical 

solution 

e T = 0.01 
(a) 

e T = 0.001 
(a) 

e T = 0.0001 
(a) 

e T = 0.00001 

0 

1.0000 

1.0000 


1.0000 


1 

2.1594 

2.1594 

2.1594 

2.1594 

2.1594 

2 

3.4522 

3.4523 

3.4522 

3.4522 

3.4522 

3 

4.8015 

4.8013 


4.8015 


4 

6.1375 

6.1372 

6.1375 

6.1375 

6.1375 

5 

7.4016 

7.0423 

7.4017 

7.4016 


6 

8.5479 

8.5481 

8.5479 

8.5479 

8.5479 

7 

9.5434 

9.5434 

9.5435 

9.5435 

9.5434 

8 

10.368 

10.367 

10.368 

10.368 


9 

11.011 

11.008 

11.009 

11.011 


10 

11.473 

11.483 

11.473 

11.473 

11.473 


Underlined digits are in disagreement with analytical solutions. 


Table II. Relative Total Flux and Computation Times for 
Analytical and Numerical Inversion Solution Methods 


[x — 100 g/cm 2 ] 


Method 

Total flux 

(relative to incident flux) 

Computation times, 
CPU sec/point 

Analytical 

11.473 

0.05 

Numerical inversion: 



£ L = 0.1 

11.500 

19.7 

£ L = 0.01 

11.480 

26.9 

£ L = 0.001 

11.474 

33.4 

e L = 0.0001 

11.473 

39.2 
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Table III. Total Flux Computations Obtained From Analytical and 
Discretized Solution Methods as Function of Atmosphere Thickness 


Thickness. 

g/cm 2 

Total flux (relative to incident flux) from — 

Analytical 

solution 

Discretized solution for number of intervals of — 

10 

20 

40 

80 

160 

20 

3.4522 

3.4688 

3.4695 

3.4634 

3.4585 

3.4556 

40 

6.1375 

5.9588 

6.0550 

6.0983 

6.1186 

6.1282 

60 

8.5479 

8.1020 

8.3219 

8.4344 

8.4911 

8.5195 

80 

10.368 

9.6984 

10.019 

10.019 

10.278 

10.322 

100 

11.473 

10.686 

11.059 

11.260 

11.365 

11.418 


Table IV. Relative Flux Differences Between Analytical and 
Multiple Collision Solutions as Function of Target Thickness 


Thickness, 

g/cm 2 

Percent 

difference 

10 

4.5 

20 

3.0 

30 

2.0 

40 

5.5 

50 

5.7 

60 

8.8 
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Table V. Values of (f>\ Determined From Analytical and Polynomial Expansion Solutions 

[J = 25; Eq = 5 GeV] 


Energy, 
MeV /nucleon 

<j > i from 
analytical 
solution 


01 from polynomial expansion solution for 

— 

One interval (single precision) with 
an allowed truncation error of — 

One interval 
(double precision) 

Four intervals 
(double precision) 

0.01 

0.001 

0.0001 

10.000 

1.3560E-03 

1.3529E-03 

1.3527E-03 

1.3527E-03 

1.3548E-03 

1.3560E-03 

18.616 

2.2176E-03 

2.2140E-03 

2.2137E-03 

2.2137E-03 

2.2176E-03 

2.2176E-03 

34.657 

3.6333E-03 

3.6252E-03 

3.6248E-03 

3.6247E-03 

3.6302E-03 

3.6333E-03 

64.520 

5.9866E-03 

5.9852E-03 

5.9846E-03 

5.9845E-03 

5.9906E-03 

5.9866E-03 

120.11 

1.0032E-02 

1.0032E-02 

1.0033E-02 

1.0033E-02 

1.0029E-02 

1.0032E-02 

223.61 

1.7661E-02 

1.7662E-02 

1.7659E-02 

1.7659E-02 

1.7656E-02 

1.7661E-02 

416.28 

3.5476E-02 

3.5471E-02 

3.5476E-02 

3.5477E-02 

3.5474E-02 

3.5476E-02 

774.96 

9.4341E-02 

9.4327E-02 

9.4339E-02 

9.4340E-02 

9.4341E-02 

9.4341E-02 

1442.7 

3.4442E-01 

3.4438E-01 

3.4442E-01 

3.4442E-01 

3.4442E-01 

3.4442E-01 

2685.8 

8.8492E-01 

8.8478E-01 

8.8491E-01 

8.8492E-01 

8.8492E-01 

8.8492E-01 
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Table VI. Dose Fraction Versus Ion Charge Number and Relative 
Integration Error for Manganese Ions Incident Upon Atmosphere 

[E 0 = 1000 MeV/nucleon] 


Ion charge 
number 

Dose fraction for relative integration error of— 

10 -2 

10“ 4 

10 -6 

10 -8 

lO -10 

10-12 

i 

5.3599E-01 

4.4442E-01 

5.3621E-01 

5.3645E-01 

5.3619E-01 

5.3452E-01 

2 

1.1879E-01 

1.3416E-01 

1.1875E-01 

1.1881E-01 

1.1876E-01 

1.1906E-01 

3 

5.9635E-02 

6.9760E-02 

5.9609E-02 

5.9608E-02 

5.9611E-02 

5.9802E-02 

4 

4.1417E-02 

4.9225E-02 

4.1398E-02 

4.1387E-02 

4.1400E-02 

4.1545E-02 

5 

3.0832E-02 

3.7123E-02 

3.0816E-02 

3.0801E-02 

3.0818E-02 

3.0933E-02 

6 

2.5981E-02 

3.1531E-02 

2.5967E-02 

2.5950E-02 

2.5968E-02 

2.6069E-02 

7 

2.0793E-02 

2.5501E-02 

2.0782E-02 

2.0764E-02 

2.0783E-02 

2.0868E-02 

8 

1.7125E-02 

2.1157E-02 

1.7115E-02 

1.7098E-02 

1.7116E-02 

1.7188E-02 

9 

1.3704E-02 

1.7033E-02 

1.3696E-02 

1.3680E-02 

1.3697E-02 

1.3756E-02 

10 

1.2350E-02 

1.5383E-02 

1.2343E-02 

1.2327E-02 

1.2343E-02 

1.2397E-02 

11 

1.0301E-02 

1.2872E-02 

1.0295E-02 

1.0281E-02 

1.0296E-02 

1.0341E-02 

12 

9.4449E-03 

1.1814E-02 

9.4393E-03 

9.4260E-03 

9.4397E-03 

9.4815E-03 

13 

8.1070E-03 

1.0156E-02 

8.1023E-03 

8.0903E-03 

8.1026E-03 

8.1386E-03 

14 

7.5260E-03 

9.4332E-03 

7.5215E-03 

7.5102E-03 

7.5218E-03 

7.5553E-03 

15 

6.5966E-03 

8.2738E-03 

6.5927E-03 

6.5826E-03 

6.5930E-03 

6.6224E-03 

16 

6.1814E-03 

7.754 7E-03 

6.1778E-03 

6.1682E-03 

6.1781E-03 

6.2056E-03 

17 

5.5054E-03 

6.9083E-03 

5.5021E-03 

5.4935E-03 

5.5024E-03 

5.5269E-03 

18 

4.7170E-03 

5.9196E-03 

4.7142E-03 

4.7068E-03 

4.7144E-03 

4.7355E-03 

19 

4.6828E-03 

5.8767E-03 

4.6800E-03 

4.6727E-03 

4.6803E-03 

4.7011E-03 

20 

4.4462E-03 

5.5800E-03 

4.4436E-03 

4.4366E-03 

4.4438E-03 

4.4636E-03 

21 

3.8863E-03 

4.8774E-03 

3.8840E-03 

3.8779E-03 

3.8842E-03 

3.9015E-03 

22 

3.5691E-03 

4.4793E-03 

3.5670E-03 

3.5613E-03 

3.5671E-03 

3.5830E-03 

23 

3.2931E-03 

4.1330E-03 

3.2912E-03 

3.2860E-03 

3.2913E-03 

3.3060E-03 

24 

3.1581E-03 

3.9635E-03 

3.1563E-03 

3.1513E-03 

3.1564E-03 

3.1705E-03 

25 

4.1970E-02 

5.2674E-02 

4.1946E-02 

4.1879E-02 

4.1947E-02 

4.2134E-02 
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Table VII. Dose Fraction Versus Ion Charge Number and Relative 
Integration Error for Neon Ions Incident Upon Atmosphere 

[E 0 — 1000 MeV/nucleon] 


Dose fraction for relative integration error of — 


Ion charge 
number 

10~ 2 

10 -4 

10~ 6 

10~ 8 

i 

5.6956E-01 

5.6942E-01 

5.7009E-01 

5.7009E-01 

2 

1.2701E-01 

1.2707E-01 

1.2687E-01 

1.2687E-01 

3 

6.3886E-02 

6.3906E-02 

6.3806E-02 

6.3806E-02 

4 

4.4419E-02 

4.4431E-02 

4.4362E-02 

4.4362E-02 

5 

3.3100E-02 

3.3108E-02 

3.3056E-02 

3.3056E-02 

6 

2.7910E-02 

2.7916E-02 

2.7873E-03 

2.7873E-02 

7 

2.2359E-02 

2.2364E-02 

2.2329E-02 

2.2329E-02 

8 

1.8430E-02 

1.8434E-02 

1.8406E-02 

1.8406E-02 

9 

1.4762E-02 

1.4767E-02 

1.4744E-02 

1.4744E-02 

10 

7.8566E-02 

7.8585E-02 

7.8463E-02 

7.8463E-02 
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Figure 10. Flux variation from analytical solution with solar source 



E, MeV/amu 


E, MeV/amu 


(a) Constant. (b) Linear. 



E, MeV/amu 
(c) Quadratic. 

Figure 11. Flux variation from power series method with solar source. 15 intervals; 2 points/interval. 



E, MeV/amu 


(b) a 0 = 0.001. 

Figure 12. Benchmark and interval calculations for <j > i 












f c (E), relative units 


10 10 ' 
E, MeV/amu 


Figure 16. Energy-dependent cross-section shape. 




E, MeV/amu 


E, MeV/amu 


(a) With energy dependence. (b) Without energy dependence. 

Figure 17. Flux variation with and without energy dependence of absorption cross section 
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